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1. Integration of BBE and relationship to equations 1-4 

1 Physical details and parameters with regard to BBE (Eq. (52)) are given elsewhere (ICRU 49). The 
critical behavior of BBE is the nonrelativistic limit, i.e. v = (or E = 0), since it produces severe 
singularities, and a cutoff at E = 1 MeV is usually assumed in Monte Carlo codes. Therefore we only 
consider here the nonrelativistic case: 

-dEldz = K-v- 2 ■ [ln(2m • v 2 / E l ) + a sheU + a Barkas + « v 2 ] 

K = {Zpl A N )-%n-q 2 -e 4 12m 

1 5 In order to avoid confusion, we note that m refers to the electron rest mass and M to the proton rest 
mass. Eq. (LI) can be integrated, if the following substitutions are carried out: 

{v 2 = 2 • £ / M; /? = 4 • m l{E L • M); E=(3 A ■ exp(-w / 2)} (L2) 

The application of Eq. (L2) to the logarithmic term of Eq. (LI) implies severe singularities, whereas 
the simultaneous integration of all terms of Eq. (LI) leads to mutual compensation of these critical 
20 terms and a cutoff is superfluous. A complete integration of Eq. (LI) with the help of Eq. (L2) and 
including relativistic correction terms (Bloch) yields the power expansions according to Eq. (1). With 
regard to water we obtain formula 2 (see also Tables LI and L2). A comparison between ICRU49 data 
and formula 2 is shown in Fig. LI; the average standard deviation amounts to 0.13 %. 



Table LI : Parameter values for Eq. (1) if E is in MeV, Ei in eV and R C sda in cm 
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Table L2: Parameter values for Eq. (2) if E is in MeV, E : in eV and R C sda in cm 
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Fig. LI. Comparison of Rcsda between formula 2 and data in ICRU49. 

5 A further possible way of integrating BBE leads to Gompertz-functions. A very low order 
approximation with the same accuracy as Eq. (2) is given by: 

N 

R CSDA = 6.8469- 1(T 4 -E + -b k exp(-^ •£„))] (L3) 

The parameters bk (dimensionless) and gk of this equation with N = 2 are (Ulmer 2007); bi = 
15.14450027, b 2 = 29.84400076, gi = 0.001260021 MeV" 1 , g 2 = 0.003260031 MeV" 1 . Eq. (L3) is also 
1 useful, since it can be inverted, i.e. E = E (Rcsda), whereas Eq. (2) leads to an ill-posed problem. The 
result is given by Eq. (3) and Table L3, if the restriction N = 5 to therapeutic protons is accounted for. 

A through and back calculation of R C sda(Eo) according to formula 2 in steps of 1 MeV up to 300 MeV 
and Eq(Rcsda) according to Eq. (3) gave a mean standard deviation of 0.1 1 %. Eq. (L3) also provides 
to use the actual (residual) energy E(z) by the substitutions Rcsda —* R-csda - z and E — > E(z). 
1 5 Equations (3 - 8) are straightforward. 



Table L3: Parameter values for the inversion Eq. (3) with N = 5, if E is in MeV, E: in eV and Rcsda in cm 
(dimension of A k : MeV/cm, p k : cm). Note: Eq. (3) only requires Ai,...,As and Pi,..,Ps. 
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Eq. (3) can be extended from water to a medium with arbitrary density p, nuclear charge Z and mass 
20 number A N by the substitutions: 



a\ =A k -(lmoyz-p-vs.i/E,)* /(A N - p j 



P k =P k -{\0-p w l\S)-{15ME i y< -A N /(p-Z) 



\(L4) 
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E(z) = (R CSDA -z)- > A\-e W [-(R CSDA -z)l /3\] 




k=\ 



Note: Ei of water amounts to 75.1 eV and p water to 1 g/cm 3 ; with respect to A N water and Z water the Bragg 
rule has been applied, i.e., A N>water /Z water = 18/10. Together with formula (L4), Eq. (3) incorporates the 
inversion of Eq. (1). With regard to heterogeneous media E(z) and dE(z)/dz can readily be determined 
in a stepwise manner, if A N , p, and Z are known (details: Ulmer and Matsinos 2010). 

The parameters of Eq. (7) and Table 1 solely refer to water. They have to be modified in similar 
fashion as by Eq. (L4): 



The powers q k (k = 1,..4) are identical with those of Table L3. 

2. Integration of the nuclear cross-section and decrease of the primary proton fluence 

In the LR of the publication (Ulmer 2007) we have given some information about a calculation model 
of the total nuclear cross-section proton - oxygen resulting from extended nuclear shell theory; a 
detailed elaboration can be found in a review (Ulmer and Matsinos 2010). The model is in very good 
agreement with measurement data of Chadwick et al 1996; the measurement data presented by Seltzer 
et al show some deviations between 30 MeV and 140 MeV. 

Fig. L2 shows that for E < 7 MeV only elastic scattering processes at the oxygen nucleus occur, and E 
= 7 MeV is the threshold energy E Th due to the Coulomb barrier. The total nuclear cross-section has a 
resonance maximum at E res = 20.12 MeV and thereafter it decreases exponentially. In the domain E > 
150 MeV and E < 270 MeV the asymptotic behavior is reached (the production of n -mesons requires 
proton energies E » 270 MeV). Fig. L2 is also very closely related to the decrease of the primary 
proton fluence. 



c p ^c p -(18/10) Z-p -{15.11 E,)^ l{A N -p w ) (k = \,..,4) 
V=>V -(10- Pw l\%)-{15.\l Erf ■ A N l{p-Z) 



(L5) 



s = 0.0922 



4 



600 
= 500 

O 

1 

5 400 
\ 300 

5 
Oi 

| 200 
| 100 







Theoretical calculation: generalized nuclear shell theory 

* Data of Berger 

• Data of Chadwick 








I # . »— ~4 , 1| 












i i i i i i i i 



30 



60 



90 



120 150 
Energy MeV 



180 



210 



240 



270 



Fig. L2. Total nuclear cross-section of the proton - nucleus (oxygen) interaction. 

By numerical evaluation, Boon (1998) presented the finding that in dependence of the primary proton 
energy the proton fluence decreases linearly. One would expect a jump, when a nuclear reaction is 
5 impossible, but due to statistical fluctuations this 'jump' is rather rounded, when the proton has lost 
most of its initial energy. An analytical integration (Ulmer 2007) also gave a linear correspondence, 
and the slope depends on E . There have been considered many attempts to fit solely this fluence 
decrease in dependence of R cst j a (see e.g. Bortfeld 1997 and references therein). Since the dependence 
of the difference E - E Th has not been accounted for in these fits, they have been made artificially 
10 complicated and are not accurate. Eq. (4) is the result of an analytical integration with & as an 
arbitrary initial fluence. Fig. L3 summarizes the results of the previous analysis (Ulmer 2007); 0o is 
normalized to 1 . For E < MeV the proton fluence again is constant without linear decrease. However, 
the range of 7 MeV protons is less than 1 mm and therefore only the roundness due to fluctuations can 
be verified. 

15 In further communications (Ulmer and Schaffner 2006, Ulmer and Matsinos 2010) we have presented 
total nuclear cross-sections for materials different from the oxygen nucleus. They have some 
importance with regard to collimator scatter and passage of protons through bones. The fluence 
decrease of primary protons according to Figures L4 - L5 is analogous to that of water (Fig. L3). A 
comparison of Figure L5 with Figure L3 shows that the slope of the fluence decrease is steeper, in 

20 particular for proton energies < 100 MeV. For copper the threshold energy E Th now amounts to 8.24 
MeV and for calcium 7.72 MeV, respectively. Due to the rather significant fluctuations we cannot 
verify a constant fluence at the end of the track. Eq. (8) can be applied to Fig. L5 with E Th = 8.24 MeV 
(copper) and E Th = 7.72 MeV (calcium). The power f for the computation of uq in Eq. (5) has to be 
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modified (f = 0.755 (copper) and f = 0.86 (calcium)): 

f(Z, A N )=a ■ A N ~' + b ■ A N ~ 2!1 + c ■ A N ~ ln + d ■ Z ■ A N ~" 2 (L7) 

Formula (L7) is closely related to the nuclear collective model (Ulmer and Matsinos 2010), where the 
5 parameters a, b, c and d are interpreted in terms of different contributions to the total cross section. 
These parameters have been determined by this model: a = -0.087660001, b = -6.379250217 , c = 
5.401490050and d = - 0.054279999. 




Formulas (8, 19) can be applied with regard to the modifications obtained by equations (1, L6 - L7). 
There are two applications, where Fig. (L5) is relevant: 1. Passage of protons through collimators. 2. 
Passage of protons through bone/metallic implants. In case 2, only a small path length has to be 
1 5 corrected. However, Figure L5 shows that the fluence decrease has also to be corrected according to 
the boundary conditions. Figures L3 and L5 do not yet provide information about the contributions 
S sp>n and S sp r . The case of nonreaction protons (nuclear potential/resonance scatter) has already been 
treated. According to Fig. 1 the contribution of reaction protons is particular important for E > 150 
MeV with increasing energy. 
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Fig. L4. Total nuclear cross-section of the proton interaction with the nuclei of some materials of interest (data 
from the Scientific Los Alamos library). 
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Fig. L5. Decrease of the fluence of 100, 200 and 250 MeV primary protons in brass (solid lines) and in Calcium 
(dashes). The Rcsda ranges are stated by perpendicular lines. 

We now present the calculation formula for this case (detailed information is given in the review of 
10 Ulmer and Matsinos 2010). Thus S sp , r is proportional to ® -2-i)-Cheavy anc ^ a function F, depending on 
some further parameters. It should be mentioned that the parameters of Eq. (L8) are not restricted to 
the oxygen nucleus. From Figure L4 and Eq. (1) corresponding parameters of some further nuclei can 
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be verified, e.g., Rc S d a , E Th and E res . We use the following definitions: 



Z R RcSDA ^ n 



r s =0.55411; r f = t, -0.000585437- (E-E res ) 
a vgl = z/^T s 2 +T m 2 +(R CSDA /4x) 2 



arg 2 = (R CSDA - z - 4i ■ n ■ z shlfl ) / ]t 2 + (r stmggle 1 1 W) 2 + {R CSDA I S ■ Inf 
a rg = ( R csDA ~ °- 5 • z shift ■ 4n-z) l{Jn ■ z Mfl ) 
The result is the following connection: 



(L8) 



S sp , r = (^o / N abs ) • O o • [2 • v ■ C heavv ■ F r +G] 



■ Q tm as {medium) ■ A N in l{Q' ot as {oxygen ) • A 



1/3 



oxygen 



) 



(p = 0.5 • (1 + erf {arg)); cp x = er/(arg 1); <p 2 = er/(arg 2) 

e Je~ X ~exp{-{z-z R ) 2 lz R 2 )l{e~ l "1) Of z^z R ) 
v l {else) y 

G = (cj • • / R CSDA ) • exp ( -(-^ 7? CSD ^ - z) 2 / z x 2 ) 



(L9) 



Eq. (L9) requires the determination of the parameters given in Eq. (L8). All parameters are stated in 
the paper, excepted E res . Thus for oxygen E res amounts to 20.12 MeV, but E res of some other materials 
(C, Ca, Cu, Zn) can be taken from Figure L4 or from a related calculation formula. With regard to the 
asymptotic cross-sections Q tot as in Eq. (L9) we may either use Figure L4 or formulas, which may be 
10 taken in a detailed discussion of equations (L8 - L9) being published (Ulmer and Matsinos 2009, 
Ulmer and Matsinos 2010). 

3. The importance of Landau-Vavilov distributions in buildup of protons and the role of 
secondary (reaction) protons 

It is a noteworthy result (Ulmer 2007) that a Boltzmann distribution function (canonical ensemble) in 
1 5 connection with a Schrodinger equation (i.e. a nonrelativistic approach) for the energy transfer from 
protons to environmental electrons yields a Gaussian convolution. This connection can briefly be 
shown as follows: 

Let O be a source function and (p an image function resulting from an exchange of particles and 
energy E ex and obeying a statistical distribution, then an exchange Hamiltonian H mutually couples the 
20 source field <1> (e.g. proton fluence) and the image field (p by the operator F H : 
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cp=F„ -0 = exp(-i/ / ^ ex ) • O 



H = -^-d 2 I dz' 

Zm 

s 2 = 2ft 2 I m • E„„ 



(L10) 



F H is an operator representation of a canonical ensemble, and the Gaussian convolution kernel K is 
obtained as the Green's function by using the spectral theorem of operators: 



F H -exp(-z- k-z)/ V^t =exp(-£ 2 ■ s 2 /4)/ 42k =^> 

CO 

K = y^- j*exp(— k 2 • s 2 /4)exp(7 • k • w)exp(— i • k • z)dk = 

-co 

K = -^exp[-((u-z)/s) 2 ] 



(Lll) 



With the help K according to Eq. (Lll) the connection between source field O and image <p can be 
represented in the more familiar fashion: 

*»•"-/*<■•-*>•*<">* ,L12) 



The energy exchange of a proton with the environmental electrons corresponds to the transfer energy 
Etransfer- Only in the thermodynamical equilibrium, which results from further collisions of electrons, 

10 we can replace E tranfer by k B -T (k B : Boltzmann constant). This aspect is also the starting point of the 
calorimetric dosimetry. Equations (L10 - Lll) directly lead to the density matrix formalism (von 
Neumann) and to an operator calculus of the quantum-statistical extension of the Feynman propagator 
theory. Instead of the spectral distribution resulting from the restriction to a simple plane wave 
representation exp(-i-k-z) it is also possible to make use of a complete power (Fourier) expansion 

15 (Ulmer2007), i.e.: 



20 



/ 2n 



exp(— i - k-z): 



I 2k 



f (k)oxp(- i - k- z) 



/(*) = 



OO 

Z 

n'=0 



a. 



k n ' 



(L13) 



This extension provides more complex spectral distributions, if F H acts on the system of functions 
provided by Eq. (LI 3) and the integration over k according to Eq. (Lll) has been carried out. It should 
also be noted that Eq. (Lll) provides general forms of the density matrix and Feynman propagators 
and may be applied to the calculation of the stopping power (see also section 2.1). The transition to 
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statistical mechanics via Green's function calculus yields an access to a Boltzmann equation and a 
connection to a paper of Sandison and Chvetsov (2000). 

If the energy transfer E transfer is repeated n-times along the proton track (i.e. difference of the actual (or 
residual) proton energy E(n) - E(n-l)), then we have to substitute s 2 — > n-s 2 . This implies that E transfer is 
5 constant along the track. This is, however, not true in the environment of the Bragg peak, and s is itself 
z-dependent, i.e. s = s(z) (see also Eq. (7), where s(z) is replaced by T strag gie). It is evident that this way 
to account for statistical fluctuations of the local energy transfer E transfer is only valid in a 
nonrelativistic approach. Thus for protons with energy » 50 MeV relativistic corrections are required, 
and only in the Bragg peak domain the significance of these corrections breaks down. 

10 If we pass from the Schrodinger equation and Boltzmann distribution function to the relativistic Dirac 
Hamiltonian and the distribution function of Fermi-Dirac statistics, i.e., a, P: Dirac matrices, p: 
relativistic momentum: 



15 



20 



H D = a p + p - mc 2 



exp((H Z5 -E jF )/E 



(L14) 



E F : energy of the Fermi edge, for the description of the fluctuations, we obtain according to Ulmer 
(2007) with the help of rather complicated calculations the convolution kernel K F and energy 
distribution function S E . They have the shape: 



K F = N f 



oo 

1=0 



B l (n, mc 2 )-H I ((u-z- z shtft (/)) / a n ) ■ exp(-(« - z - z shtft (I)) 2 1 2a n ) 



(L15) 



S E = N, ■ cxp(-(E n (k) - E Averamn ) 2 /2a E (n) n )) • 

CO 



1=0 



(L16) 



77/ 2 4 , 4-2/2 2 

t = m c +n k c 



N f is a normalization factor and the repetition factor n has the same meaning as above. Note that 
E A verage,n results from the energy of the Fermi edge E F , and by restricting Eq. (LI 6) to the lowest order, 
i.e., / = and all higher order terms are omitted, Bohr's classical formula of energy fluctuations is 
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obtained. 

Now E Av era g e and z s hift depend in every order on the energy of the Fermi edge E F . The expansion 
coefficients bi and Bi are determined by a uniformly convergent expansion of the sech-function in 
terms of a Gaussian and Hermite polynomials Hi. However, for the purpose of therapeutic protons we 
do not need fully exploit these expansions. In a low order approximation, equations (LI 5- LI 6) imply 
a Vavilov distribution, containing Landau tails. For very low energies (nonrelativistic limit: E — > 0), 
the limit of a Gaussian kernel is obtained again. In particular, a low-level approximation of Eq. (LI 6) 
is of interest for analytical convolutions ('generalized Gaussian kernel') to account for deviations of 
symmetrical fluctuations. A glance to Fig. L6 shows that E max of the energy transfer of 250 MeV 
protons amounts to ca. 617 keV and decreases along the proton track. The straight-line represents the 
nonrelativistic determination of E max . The contributions I Lan i and I Lan2 in section 2.6 represent 
integrations of terms up to the order 2 of Eq. (LI 5). 
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Fig. L6. Maximum transfer energy E max of protons in dependence of the energy (dots: ICRU49). The straight line 
below indicates a linear connection between proton energy and maximal energy transfer resulting from a 
nonrelativistic calculation. 



In order to obtain a fast approximation up to order 2 of Landau tails according to equations (28 - 29), 
we have assumed that in the environment of the impinging surface the first and second order terms are 
leading, whereas in the domain of decreased proton energies only terms up to order 1 are significant 
(they are certainly also of interest at surface). All necessary calculation parameters have been 
determined by fits to the more general theory as outlined above and by fits to Monte Carlo (GEANT4), 
since the code has also the option to use a model of Vavilov distributions. 
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In the following, we abbreviate T Lan i with Xi and T La n2 with N fl and Ne are suitable proportionality 
factors . The determination of formula 28 results from the first order Hermite polynomial. In principle, 
we have to evaluate the convolution integral 



^1 =H fl ■T;-— R J (u-z^-expi-iu-z-z^f/T*] (L17) 



5 This integral can be evaluated by the substitution: 
u' = (u-z-z Ls )/t 1 (LI 8) 

I i m \= N f x -^T' \(u'-T l +z)-QM-u a )du ' (L19) 



■\ ■ R c S da J 



-(z is +z)/r, 



Eq. (L1 9) simply yields error functions and Gaussians, which we write in the form l Lan i =T er f + T gauss : 
T erf= N fi ■^■[erf(z+z Is )/T l )+erf((R csda -z-z Is )/r l )] (L20) 

10 T 8ma =N fi •z l /(^r-R csda )[ Q M-(z + z L f I T x 2 ))- QW {-{R csda -z-z L f I r 2 ))] (L21) 
The Hermite polynomial of second degree (unnormalized) reads: 
H 2 (£) = 4-{ 2 -2 (L22) 

We have simply to evaluate 



R 

15 ^„2-4.7V /2 .-^.i. J ?? ia „ 2 - 2 . I (u-z^) 2 -esp(-(u-z-z u y I t 2 z )du + T c (L23) 



J' 





The term T c refers to the convolution of the factor '- 2' in H 2 according to relation L22, which is only a 
'standard problem'. We replace Ti by x 2 in the substitution LI 8: this yields the following evaluation 
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\ t-- \ 2 , 2 , _2 



(w 1 -r 2 +z + 2 ■ w '-z) ■ exp(— w 1 )du'+T c 

(z+z l J/t 2 



T c =-2-N f ■ + 

C J2 4ti 



exp(-w ' 

-(z+Zi,)/r 2 



(L24) 



Thus the terms containing z , 2-u'-z and T c have already been solved. The remaining term containing 
u' 2 -x 2 2 can be integrated with the help of the chain rule, i.e.: 



b 



u a ■ exp(-u' 2 )du' =-0.5 -(b-exp(-b 2 )-a- exp(-a 2 )) + T erf (a,b) (L25) 



5 The integration boundaries a, b are stated in Eq. (L26). 

T erf = 0.5 • f ■ [erf (b) -erf (a)] j ^ 

« = ( z + z av)/^2 ; b = (R Lan2 -z-z u )It 2 \ 

It seems noteworthy to add some comments referring to equations (L17 - L26). In equations (28 - 29) 
we have stated not all terms as required by the above evaluation. The reason is that we have carried 
out a selection and accounted only for those terms, which lead to significant contributions. Terms of 
10 the order < 0.2 % have been neglected. This is justified, since the overall contributions of Landau tails 
are small in the energy domain of therapeutic protons and saving computation time should also be 
accounted for. 



Carlsson et al (1977) have reported a buildup of the Bragg curve of 185 MeV protons (E max of the 
1 5 incident beam: ca. 446 keV) and interpreted this effect as a sole result of secondary protons, since 
these protons are rather generated along the proton track than at surface. Following aspects should be 
noted: 

1 . We have already verified that secondary reaction protons significantly contribute to the buildup of 
monoenergetic protons. With regard to polychromatic protons the spectral parameter x in implies an 

20 enhancement of the buildup by accounting for the Landau tail. Further modulations are possible by 
collimator scattering occurring in broad proton beams and depending on the field-size. 

2. The question also arises, whether the transport of 8-electrons is responsible for the buildup similar 



13 



to the Compton scatter of photons. However, E^ of 250 MeV protons amounts to 617 keV (Fig. L6) 
and with regard to the cited 185 MeV protons this energy is lower. The range of these electrons is too 
small to explain a buildup. In order to produce E max of the order 4 MeV, the proton energy should 
amount to ca. 1 GeV. 

5 3. A possible contribution of the y-quanta resulting from P+-decays of heavy recoils and y-quanta by 
the annihilation of positrons cannot be significant, since the interaction processes are rare enough. 

In agreement with GEANT4 the depth dose curve of primary protons (250 MeV) shows a valley in the 
middle part of the plateau resulting from the corresponding fluence decrease (Fig. 2). If the transport 
of all secondary protons is included, the total depth dose curve does not show this valley. A 

1 comparison with the results of Medin et al (1997) is noteworthy. If the transport of secondary protons 
is only partially accounted for or even omitted (PTRAN), this valley can also be recognized at the total 
depth dose curve. In addition, it is necessary to include the 'secondaries' in an accurate manner. 
According to the cited authors PTRAN leads for 200 MeV protons to a dose contribution of ca. 10 % 
at the depth z = 20 cm, whereas very early theoretical calculations of Zerby et al (1965) provided 17 % 

1 5 at the depth z = 20 cm. This is in agreement with Fig. L3 and may be calculated with the help of 
section 2.5. In order to be consistent with the total nuclear cross-section (Fig. L2), the resulting Fig. L3 
and the classifications in section 2.5, we consider all protons as secondary protons, which result from 
the total nuclear cross-section. The so-called 'reaction protons (sp,r)\ which stand in a close 
relationship to the heavy recoils, are very dominant for E > 150 MeV with increasing tendency. Their 

20 depth dose curve certainly shows a maximum along their track, but not a typical Bragg peak due to the 
broad spectral distribution (Fig. 1). The secondary (nonreaction) protons lead to a Bragg peak shifted 
slightly to a lower energy due to further losses of energy (resonance scatter). In various publications 
on therapeutic protons this distinction has not clearly been pointed out. 

In relationship to figures 1 - 4, which deal with the role of reaction protons and the Landau tail of 250 
25 MeV monoenergetic protons, we present some further calculations of primary protons 
(monoenergetic) and comparisons with GEANT4 calculations (Fig. L7). If we glance at this figure, 
which shows complete Bragg curves, we are hardly able to verify significant differences between 
GEANT4 and theoretical calculations. Therefore, we additionally present the buildup domain (see Fig. 
L8) and the behavior of the Bragg curves in the region with z > 15 cm (Fig. L9). Thus Fig. L8 clearly 
30 shows the small influence of the asymmetric energy transfer (Landau tail) to the buildup in the high 
energy domain. Figures L9 and L10 yield an indication of the role of cutoffs in GEANT4 due to the 
severe singularities in the environment of the distal end. As already pointed out (section LR.l and 
Ulmer 2007) these singularities emerge from the behavior of BBE in the low energy domain, i.e., if 
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either E — > or if the logarithmic term becomes 0. By that, the singularity problem is a general feature 
of Monte-Carlo codes. Resulting from the energy/range straggling (T strag gie) and the cutoff, the peak 
height calculated by GEANT4 is slightly lower than by the presented theoretical calculation. In all 
other remaining domains the difference between GEANT4 and theoretical calculation is extremely 
small, i.e., less than 0.2 %. With reference to Figure L8 it should be mentioned that the small buildup 
effect due to the Landau tail disappeared in GEANT4, if the Vavilov distributions was reduced to a 
Gaussian one as mentioned in section 2.6. 
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Fig. L7. Comparison of calculated stopping powers of monoenergetic primary protons with GEANT4 (dots). 
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Fig. L9. Section of Fig. L7 in order to indicate the differences at the distal end. 
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Fig. L10. 150 MeV protons at the distal end (solid: calculation, dashes: GEANT4). 



4. Multiple scatter procedure of primary and secondary protons 

5 In the present implementation of Eclipse, the lateral scatter of protons is treated by an approximate 
version of the multiple scatter theory (Bethe 1953, Moliere 1955, Gottschalk et al 1993, Gottschalk et 
al 1999). However, we improved the Highland approximation (Highland 1995 and 1997) by the use of 
two Gaussian kernels in order to describe the lateral tail of primary protons and secondary 
(nonreaction) protons (sp,n) more accurate. GEANT4 indicated that two Gaussian kernels are 

1 sufficient for primary protons. Thus the first Gaussian accounts for the inner part of multiple Moliere 
scatter, which is steeper than the Highland approximation, whereas the second Gaussian has a much 
larger half-width to approximately describe the tail. The Highland approximation assumes a slightly 
broader half-width than necessary for the inner part in order to account partially for the tail. For 
secondary (reaction) protons (sp,r), we restricted the lateral kernel to one modified Gaussian due to the 

1 5 significantly smaller contribution of them. Thus, the calculation model of the lateral kernel for water is 
given by: 

r 2 r 2 

K lat , prim (r,z) = -^-e~^ + e'^f (L27) 
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We use a weight coefficient C of 0.96 for the contribution of the main Gaussian. The calculation of 
Ti at (z) (inner part) and Ti atLA (z) (large angle) is carried out as follows: 

T lat {z)=f- 0.626 - W e(z) 
/ = 0.9236 



1.61418-z 
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(E 7176.576)^; p = l.5 + 



[0.00150 -(l 76.576 -E ), if E < 176.576 
10.03104 • -176.576, if E >176.576 



(L28) 



The parameter x dista i will be discussed at the end of this section. A good approximation for a model 
with a single Gaussian for primary protons can be obtained from the above equations by substituting 
C and f by 1. In agreement with results of Gottschalk et al 1993 the change from water to other 
materials is obtained by scaling of Q(z) with the help of a different value for Rcsda- The error function 



18 



term models the Gaussian distribution of the stopping distribution due to range straggling. The protons 
which have undergone only very small angle scattering are closer to the central axis of the beamlet and 
travel slightly further. 

A fit of Monte Carlo results shows that Ti at ,LA(z) can be determined by the kernel: 
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(L29) 



The scaling properties of Q(z) still hold, since only the ratio z/R C sda enters Eq. (L29). For the 
secondary reaction protons and heavy recoils, we use separate Gaussian kernels, i.e.: 



(L30) 



^heavy( Z )^ T heavy( E ( Z ))\ 

Q(z) and are the same as defined in Eq. (L28), but T he av y (E(z)) is given by the same Eq. (24) as 
above except that it depends on the local energy E(z) instead of T he avy(Eo). 



It might appear that Eq. (L27) and the calculation of x max according to Eq. (L28) are just obtained by 
fitting methods. However, the spatial behavior of the multiple scatter theory can be written by a 
Gaussian and Hermite polynomials H 2n : 
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(L31) 



N = \l{47t-lat(z)); r 2 =x 2 +y 2 
The parameters of the Hermite -polynomial expansion are: a = 0.932; a 2 = 0.041; a 4 = 0.019; a 6 = 
0.008. It is now the task to determine a linear combination of two Gaussians according to Eq. (L27) 
with different half-widths by a least square -fit (standard task in many problems). This way, we are 
able to go beyond the Highland approximation. 

The calculation of is carried out in the following way: The differential cross-section is given by: 



dlat(z) 2 _ dlat 
dz 

2 
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(L32) 



dE E 2 

In our calculations, we only used a W ater- a Me dium is proportional to Z-p/A M , i.e. nuclear charge, mass 
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density and mass number. An evaluation of E(z) and dE/dz in Eq. (L32) above can be simplified, when 
using the Bragg-Kleeman rule: 



R CSDA =0.00259 -E P 



0.00259 • e{z) } 
The inversion of Eq. (L33) is given by: 



^csda z 



(L33) 



E(z) = 
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CSDA 



— Z 



0.00259 



(L34) 



An accurate application of this rule requires consideration of p depending on E . p is in the order of 1.7 
- 1.8 according to Ulmer (2007). With the help of Eq. (L34) dE(z)/dz can be computed and the 
evaluation of Eq. (L32 ) yields: 
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(L35) 



1 One might expect that the lateral scatter functions according to equations (L28 - L29) continuously 
increase up to z = R csda . However, this assumption would be valid, if the energy spectrum for scattered 
protons would be identical at a depth z independent of the scatter angle and fluctuations due to the 
energy/range straggling T stra ggie and x in , respectively. In fact, there are small fluctuations of the lateral 
scatter functions along the proton tracks. In particular, from Bragg peak down to the distal end there is 

1 5 a significant difference between those protons, which have only undergone small angle scatter in this 
domain, and those with larger scatter angle. The latter protons have deposited their energy in an 
oblique path, and therefore they stop earlier and cannot reach the distal end of the Bragg curve. It is 
clear that the scatter functions for primary and secondary protons (x ta , Xi at>LA and x sp ) depend on x strag gi e 
and x in , which are the origin of these fluctuations and yielding the significant change of the energy 

20 spectrum at the end of the proton tracks. In order to describe this behavior by a mathematical model, 
we prefer to use a Gaussian convolution, which is certainly justified in the domain of Bragg peak (low 
energy region of the proton tracks). As an example we use the function Q(z) in Eq. (L28), which 
determines both Xi at and x sp . We denote the fluctuation parameter by Xdi sta i; the connection to straggle 
parameters will be considered thereafter. 



25 The crude model assumes: 
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fio = 



e xp(1.161418-z/7? cm ,)-l \ 
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1 (i/ z>R CSDA ) 



(L36) 



The more realistic model taking all the arguments with regard to fluctuations into account is obtained 
by: 



D 

00) Jeo(^)exp[-(z-^) 2 /^ ta /]^ (L37) 

— GO 

5 Using methods, often applied in this paper, we obtain the modified Eq. (L28) for Q(z) instead of Qo(z) 
according to Eq. (L37): 

n(-\- eX P( U61418 ' Z l R csda) "I I \\s pr f((T> - V \It Yl 

^ ' exp(l. 161418) - 1 2 L i_t ~^V U^csc/a ^ ) 1 L distal ) \ (L38) 

This result implies that Q(z) increases exponentially along the proton track, as long as the error- 
function is 1. Only in the environment of z = R csda Q(z) decreases rapidly. However, this depends 
1 actually on x. The connection between x and the above stated convolution parameters for energy/range 
straggling is a principal question. One might assume that for proton pencil beams with energy/range 
straggling we can set: 

2 2 2 

T distal =T straggle +T in (L39) 

This assumption might be reasonable, since the proton history resulting from the beam-line and 
15 expressed by x in should be accounted for. As already pointed out we have made use of GEANT4 
results for the adaptation of the scatter functions with regard to monoenergetic protons. It turned out 
that the best adaptation in the domain from the Bragg peak to the distal end can be obtained, if we put 
for monoenergetic protons: 

2 2 2 

T distal = T straesle + 1 Ranee (L40) 



' distal straggle 

20 The mean standard deviation for protons with E = 50 MeV up to 250 MeV in intervals of 25 MeV 
amounts to 2.6 %. x Range is given by Eq. (17). Therefore it may be justified to modify Eq. (L40) by: 
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2 2 



T =T . , +T. + T a (L41) 

straggle in Range v ' 

TRange (Eq. ( 1 7)) is also a convenient parameter to include the initial phase space of beamlets, since the 
radial energy dependence of a beamlet (i.e. decrease of the energy) implies varying ranges relative to 
25 the central ray. With regard to the evaluation of all scanning data of an Accel machine the corrections 
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by fitting x Ra n g e to account for initial phase space properties amounted to x Range ,fi t = x Range (1 + e) with 8 
< 0.12. Thus by the substitution x Range — » x Rang e,fit a slight correction of Eq. (L41) is obtained. 

An important aspect is the handling of the above formulas in the case of media different from water. 
According to findings of Gottschalk et al (1993) we are able to represent lateral scatter functions in a 
dimensionless manner approximately valid for arbitrary media by the substitutions x' = Xi at /x max and s = 
z/Rcsda, i.e., x' =x'(s) should be independent on the medium. Gottschalk et al demonstrated that this 
relationship approximately holds even for lead. Figure Lll presents the findings of Gottschalk et al 
and incorporates an essential tool with regard to scatter in heterogeneous media. An evaluation of x max 
with tools worked out in a review (Ulmer and Matsinos 2010) provided the scaling factors: 
w(calcium) = x inax (water)T.0281; x^copper) = x max (water)T.0394; x max (lead) = x max (water)T.0826. 
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Fig. Lll. Scaled representation of the lateral scattering function of water, copper and lead (Gottschalk et al 
1993) 



1 5 5. Calculation of the absolute normalization factor 

The normalization factor N abs , which appears in all contributions, provides normalization of the Bragg 
curves to the absolute dose in MeV/cm (in water). All distances are stated in units of cm and the 
energy in MeV. N abs is important with regard to the monitor unit calculation, A is the area integral of 
the complete integral over S(z), if S(z) is given by relative values (1 or 100 % at the Bragg peak). All 
20 necessary integrals and can be carried out analytically. 
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(L42)+ 



++ In order to make apparent the integration procedure we state an example by Eq. (L48 - L49) at the 
end of this section. The recoil protons are treated here as primary protons. The integral over the 
contribution of primary protons is split into two parts for computational reasons. A pp i is calculated 
under the assumption, that there is no loss of primary protons, i.e. <& pp = <J> - 
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(L43) 



A PP 2 represents the energy lost by the primary protons due to nuclear reactions: 



uq 



ppl ' 



■Jx ■ 1 + erf 



( R 


•\ 


\ 




csda 




^st 








aggie j 


) 



2 A 



+ T 



simple 



+ C 2- R csda +2-C 3 



PE-n 



PE-n- R 



J 

PE-a-R, 



•1 + e 



+ • C4 ' R CS( i a +-£"() • C !ntl \ 'R Jlin \ • C/„„9 •-/^/, 



(L44) 



A sp is obtained by an integration of the contribution of secondary protons S sp . A spn is somewhat 
1 smaller than A pp2 , since it does not include the energy carried away by photons or neutrons or used in 
nuclear reactions. 
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(L45) 
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Rcsda = R-csda ~ Z shift 
R~Lan\ ~ R-Lan\ ~ Z shift 
^Lanl ~~ ^Lan2 ~ Z shift 



(L46) 



Arp and A sp r yield a rather complicated equations; they can be simplified by application of the mean 
value theorem of the integral calculus: 

= C heavy • (R csda — V2 -T heavy —3mr Range ) 



sp,r 



2 • C heavy ■ v ■ 4n ■ (R csda - z m ■ 4n~) + \ ■ 4n ■ c x ■ z sUft I R csda 



(L47) 



f The integration procedure underlying Eq. (L42) is shown by the term C 2 of Eq. (14): 



(limz ^co) C 2 -j[\ + erf((R csda -z)/T)]dz = 

= C 2 -[2-z -t-(u 2 ■erf{u 2 )-u x ■erf(u l )) + j^-(eyp(-u 2 2 )-exp(-u l 2 ))] = 2-C 2 ■ R csda 

U 2=( R csda- Z o) /T '> U l=( R csda +Z o) /r > Z ~> 00 

In order to evaluate the integral over erf((R csda - z)/x) of formula (L48) we have to perform the 
substitution u = (R cs d a - z)/x and to use the chain rule: 



(L48) 



10 I erf{u)du=u-erf{u) + -j=-Q-m(-u 2 ) (L49) 
J ^ 

Since erf((R csda - z )/x) rapidly assumes -1 for z — > °o and erf((R csda +z )/x) — >■ 1, we finally obtain 
2-C 2 -R C sda as a contribution to A pp i; the Gaussian terms in Eq. (L48) mutually cancel for z — > oo. The 
evaluation of all other terms leading to equations (L42 - L47) follows the same principle. 



